R Notebook

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.0 
✔ readr   2.1.2      ✔ forcats 0.5.2 ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(geojsonio)

Registered S3 method overwritten by 'geojsonsf':
  method        from   
  print.geojson geojson

Attaching package: ‘geojsonio’

The following object is masked from ‘package:base’:

    pretty

library(broom)
library(sf)

Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(osmdata)

Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

library(ggplot2)
library(dplyr)
library(ggnewscale)
library(RSocrata)
library(rstudioapi)
library(leaflet)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(shiny)

Attaching package: ‘shiny’

The following object is masked from ‘package:geojsonio’:

    validate

library(bslib)

Attaching package: ‘bslib’

The following object is masked from ‘package:broom’:

    bootstrap

The following object is masked from ‘package:utils’:

    page

library(RColorBrewer)
#set working directory
current_path <- rstudioapi::getActiveDocumentContext()$path
setwd(dirname(current_path))
getwd()

[1] "/Users/xxxxoxygene/Downloads/Columbia University/Fall2022/STAT4243/Project 2/fall2022-project2-group2/doc"

#### Load the Restaurant dataset
data =  read.socrata(
  "https://data.cityofnewyork.us/resource/43nn-pn8j.json",
  app_token = "zTRehp1897SQtpYtBiIOUMfR4"
)

data$year <- format(data$inspection_date,"%Y") ## Extract Years

###filter the dataset
df =
  data %>%
  filter(data$year >= 2019 & zipcode!="" & dba != "") %>%
  mutate(grade=replace(grade, grade == "", NA))

head(df)

# Read the geojson file containing Spatial Info
spdf_file <- geojson_read("../data/zip_code_040114.geojson", what = "sp")

stats_df = spdf_file@data

# Convert it to a spatial data frame, with zip code as index
spdf_data <- tidy(  spdf_file,
  region="ZIPCODE"  # Use ZIPCODE variable as index, the index will be named "id"
)

Section II: interactive map

##### None Interactive map (Population by region)
ggplot() +
  geom_polygon(data=spdf_data %>%
                 left_join(stats_df, c("id"="ZIPCODE")),
               aes(x=long,
                   y=lat,
                   group=group,
                   fill=POPULATION),
               color="white",
               size=.2) +
  theme_void() +
  coord_map() +
  scale_fill_distiller(palette = "YlGnBu", direction = 1) +
  labs(title="Population in New York City",
       subtitle="Neighborhoods are filled by population",
       fill="Population")

#### Number of restaurant per ZIPCODE
Num_Rest_Code =
  df%>%
  group_by(zipcode, dba, latitude, longitude)%>%
  count() %>%
  group_by(zipcode)%>%
  count()

Critical_2019_by_Code = 
  df%>%
  filter(year == 2019)%>%
  group_by(zipcode)%>%
  summarize(Total = n())
  
Critical_2020_by_Code = 
  df%>%
  filter(year == 2020)%>%
  group_by(zipcode)%>%
  summarize(Total = n())
  
Critical_2021_by_Code = 
  df%>%
  filter(year == 2021)%>%
  group_by(zipcode)%>%
  summarize(Total = n())
  
Critical_2022_by_Code = 
  df%>%
  filter(year == 2022)%>%
  group_by(zipcode)%>%
  summarize(Total = n())


Critical_spdf_file_2022 = spdf_file
Critical_spdf_file_2022@data =
Critical_spdf_file_2022@data %>%
                 left_join(Critical_2022_by_Code, c("ZIPCODE"="zipcode"))



Critical_spdf_file_2019 = spdf_file
Critical_spdf_file_2019@data =
Critical_spdf_file_2019@data %>%
                 left_join(Critical_2019_by_Code, c("ZIPCODE"="zipcode"))


Critical_spdf_file_2020 = spdf_file
Critical_spdf_file_2020@data =
Critical_spdf_file_2020@data %>%
                 left_join(Critical_2020_by_Code, c("ZIPCODE"="zipcode"))

Critical_spdf_file_2021 = spdf_file
Critical_spdf_file_2021@data =
Critical_spdf_file_2021@data %>%
                 left_join(Critical_2021_by_Code, c("ZIPCODE"="zipcode"))


critical_violations = list(Critical_spdf_file_2019,Critical_spdf_file_2020,Critical_spdf_file_2021,Critical_spdf_file_2022)
names(critical_violations) <- c("2019","2020","2021","2022")



nc_pal= colorNumeric(palette="YlOrBr", domain= Critical_spdf_file_2019@data$Total, na.color = 'transparent')



leaflet()%>%
  addProviderTiles("CartoDB")%>%
  #### First Layer of PolyGons
  addPolygons(
    data = Critical_spdf_file_2022 ,
    weight = 0.5,
    color = "black",
    stroke=TRUE ,
    opacity = 1 ,
    fillColor = ~nc_pal(Total),
    label = ~paste0 ('Total Critical Violation : ' , Total),
    group = '2022',
    highlight = highlightOptions(weight  = 3, color = "red", bringToFront =  T)
    ) %>%
  
  #### Second Layer of PolyGons
    addPolygons(
    data = Critical_spdf_file_2021 ,
    weight = 0.5,
    color = "black",
    stroke=TRUE ,
    opacity = 1 ,
    fillColor = ~nc_pal(Total),
    
    label =~paste0 ('Total Critical Violation : ' , Total),
    group = '2021',
    highlight = highlightOptions(weight  = 3, color = "red", bringToFront =  T)
    ) %>%
    addLayersControl(overlayGroups = c("2022", "2021"))%>%
  
  
  
    addLegend( pal=nc_pal, values= Critical_spdf_file_2022$Total, opacity=0.9, title = "Critical", position = "bottomleft" )

Warning: Some values were outside the color scale and will be treated as NAWarning: Some values were outside the color scale and will be treated as NA

NA

Total_violation = 
  df%>%
  group_by(zipcode)%>%
  summarize(Total = n())


Total_2019_by_Code = 
  df%>%
  filter(year == 2019)%>%
  group_by(zipcode)%>%
  summarize(Total = n())
  
Total_2020_by_Code = 
  df%>%
  filter(year == 2020)%>%
  group_by(zipcode)%>%
  summarize(Total = n())
  
Total_2021_by_Code = 
  df%>%
  filter(year == 2021)%>%
  group_by(zipcode)%>%
  summarize(Total = n())
  
Total_2022_by_Code = 
  df%>%
  filter(year == 2022)%>%
  group_by(zipcode)%>%
  summarize(Total = n())





##### Join datasets
Total_spdf_file_2022 = spdf_file
Total_spdf_file_2022@data =
Total_spdf_file_2022@data %>%
                 left_join(Total_2022_by_Code, c("ZIPCODE"="zipcode"))



Total_spdf_file_2019 = spdf_file
Total_spdf_file_2019@data =
Total_spdf_file_2019@data %>%
                 left_join(Total_2019_by_Code, c("ZIPCODE"="zipcode"))


Total_spdf_file_2020 = spdf_file
Total_spdf_file_2020@data =
Total_spdf_file_2020@data %>%
                 left_join(Total_2020_by_Code, c("ZIPCODE"="zipcode"))

Total_spdf_file_2021 = spdf_file
Total_spdf_file_2021@data =
Total_spdf_file_2021@data %>%
                 left_join(Total_2021_by_Code, c("ZIPCODE"="zipcode"))


total_violation <- list(Total_spdf_file_2019,Total_spdf_file_2020,Total_spdf_file_2021,Total_spdf_file_2022)
names(total_violation) <- c("2019","2020","2021","2022")


##### colors
nc_pal= colorNumeric(palette="YlOrBr", domain= Total_spdf_file_2019@data$Total, na.color = 'transparent')


leaflet()%>%
  addProviderTiles("CartoDB")%>%
  #### First Layer of PolyGons
  addPolygons(
    data = Total_spdf_file_2022 ,
    weight = 0.5,
    color = "black",
    stroke=TRUE ,
    opacity = 1 ,
    fillColor = ~nc_pal(Total),
    label = ~paste0 ('Total Critical Violation : ' , Total),
    group = '2022',
    highlight = highlightOptions(weight  = 3, color = "red", bringToFront =  T)
    ) %>%
  
  #### Second Layer of PolyGons
    addPolygons(
    data = Total_spdf_file_2021 ,
    weight = 0.5,
    color = "black",
    stroke=TRUE ,
    opacity = 1 ,
    fillColor = ~nc_pal(Total),
    label =~paste0 ('Total  Violation : ' , Total),
    group = '2021',
    highlight = highlightOptions(weight  = 3, color = "red", bringToFront =  T)
    ) %>%
    addLayersControl(overlayGroups = c("2022", "2021"))%>%
  
  #####Third layer
    addPolygons(
    data = Total_spdf_file_2020 ,
    weight = 0.5,
    color = "black",
    stroke=TRUE ,
    opacity = 1 ,
    fillColor = ~nc_pal(Total),
    
    label =~paste0 ('Total  Violation : ' , Total),
    group = '2020',
    highlight = highlightOptions(weight  = 3, color = "red", bringToFront =  T)
    ) %>%
  
  ####Fourth
    addPolygons(
    data = Total_spdf_file_2019 ,
    weight = 0.5,
    color = "black",
    stroke=TRUE ,
    opacity = 1 ,
    fillColor = ~nc_pal(Total),
    label =~paste0 ('Total  Violation : ' , Total),
    group = '2019',
    highlight = highlightOptions(weight  = 3, color = "red", bringToFront =  T)
    ) %>%
   
  
    addLayersControl(overlayGroups = c("2022", "2021",'2020', '2019'))%>%
    addLegend( pal=nc_pal, values= Total_spdf_file_2022$Total, opacity=0.9, title = "Count of Total Violation", position = "bottomleft" )

Warning: Some values were outside the color scale and will be treated as NAWarning: Some values were outside the color scale and will be treated as NA

violations <- list(total_violation,critical_violations)
names(violations) <- c("Number of Total Violations","Number of Crital Violations")

restaurant info

for (y in c("2019","2020","2021","2022")){
  score_year[y] <- score[score$year==y,]
}

Warning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement length

Section III: Rshiny App

nc_pal= colorNumeric(palette="YlOrBr", domain= Total_spdf_file_2022@data$Total,na.color = 'transparent')
ui <- navbarPage(
 theme = bs_theme(bootswatch = "litera"),
  "Food Inspectation",
  tabPanel("Introduction"),
  tabPanel("Static Plots"),
 
  navbarMenu("Interactive Plots",
             tabPanel("Interactive Map",
           fluidRow(column(6,selectInput("type","Type of Violations:",
                                    c("Number of Total Violations",
                                      "Number of Crital Violations"))),
                    column(6,selectInput("time","Year:",
                                    c("2019","2020","2021","2022")))),
          fluidRow(leafletOutput("map",height = 1000))),
          
          
            tabPanel("Comparison between Years",
                      fluidRow(column(4,selectInput("type_comp","Type of Violations:",
                                    c("Number of Total Violations",
                                      "Number of Crital Violations"))),
                    column(4,selectInput("time1","Year:",
                                    c("2019","2020","2021","2022"))),
                    column(4,selectInput("time2","Year:",
                                    c("2019","2020","2021","2022"),selected = "2020"))),
                    fluidRow(column(6,leafletOutput("map_comp1",height=600)), column(6,leafletOutput("map_comp2",height=600))))),
 
  tabPanel("Reference")
)
server <- function(input, output,session){
  #interactive map
  output$map <- renderLeaflet({
    leaflet()%>%
  addProviderTiles("CartoDB")%>%
  #### First Layer of PolyGons
  addPolygons(
    data = violations[[input$type]][[input$time]],
    weight = 0.5,
    color = "black",
    stroke=TRUE ,
    opacity = 1 ,
    fillColor = ~nc_pal(Total),
    label = ~paste0 ('Total Violations : ' , Total),
    group = '2022',
    highlight = highlightOptions(weight  = 3, color = "red", bringToFront =  T)
    )%>%
    addLegend( pal=nc_pal, values= violations[[input$type]][[input$time]]$Total, opacity=0.9, title = "Count of Total Violation", position = "bottomleft" )
    })
  
  
  
  
#interactive map compared by year  
    output$map_comp1 <- renderLeaflet({
      leaflet()%>%
      addProviderTiles("CartoDB")%>%
      addPolygons(
      data = violations[[input$type_comp]][[input$time1]],
      weight = 0.5,
      color = "black",
      stroke=TRUE ,
      opacity = 1 ,
      fillColor = ~nc_pal(Total),
      label = ~paste0 ('Total Violations : ' , Total),
      group = '2022',
      highlight = highlightOptions(weight  = 3, color = "red", bringToFront =  T)
      )%>%
      addLegend( pal=nc_pal, values= violations[[input$type]][[input$time1]]$Total, opacity=0.9, title = "Count of Total Violation", position = "bottomleft" )
    })
    output$map_comp2 <- renderLeaflet({
      leaflet()%>%
      addProviderTiles("CartoDB")%>%
      addPolygons(
      data = violations[[input$type_comp]][[input$time2]],
      weight = 0.5,
      color = "black",
      stroke=TRUE ,
      opacity = 1 ,
      fillColor = ~nc_pal(Total),
      label = ~paste0 ('Total Violations : ' , Total),
      group = '2022',
      highlight = highlightOptions(weight  = 3, color = "red", bringToFront =  T)
      )%>%
      addLegend( pal=nc_pal, values= violations[[input$type]][[input$time2]]$Total, opacity=0.9, title = "Count of Total Violation", position = "bottomleft" )
    
    
    })
}

shinyApp(ui,server)